$ontext
This code is for the following paper:
    Cai, Yongyang, William Brock, and Anastasios Xepapadeas (2022).
    Climate Change Impact on Economic Growth: Regional Climate Policy under Cooperation and Noncooperation.
    JAERE, forthcoming

Author: Yongyang Cai, The Ohio State University
Date: August 11, 2022
$offtext

set n sets of scenarios / RCP26, RCP45, RCP6, RCP85 /;

* years in MAGICC6 data
SET  t                 Time periods                     / 1765*2310 / ;

* years in CMIP5 data
set tt(t) / 2006*2100 /;

* years for calibration
set ttt(tt) / 2015*2100 /;


$offlisting

Table Etotal(t,n) 
$ondelim
$include RCP_TotalE.csv
$offdelim

Table CA(t,n) 
$ondelim
$include RCP_MAT.csv
$offdelim

$onlisting
;
 
* 1 ppm = 2.13 GtC
CA(T,n) = CA(T,n)*2.13; 


parameters
** Carbon cycle
* Initial Conditions
    mat0   Initial Concentration in atmosphere 2015 (GtC)        /851   /
    mu0    Initial Concentration in upper strata 2015 (GtC)      /460   /
	ml0    Initial Concentration in lower strata 2015 (GtC)      /1740   /
    mateq  Equilibrium concentration atmosphere  (GtC)           /588     /
    mueq   Equilibrium concentration in upper strata (GtC)       /360    /
    mleq   Equilibrium concentration in lower strata (GtC)       /1720   /

* Flow paramaters
    b12          Carbon cycle transition matrix                     
	b11          Carbon cycle transition matrix
 	b21          Carbon cycle transition matrix
	b22          Carbon cycle transition matrix
	b23          Carbon cycle transition matrix		
 	b32          Carbon cycle transition matrix
	b33          Carbon cycle transition matrix
;

b12 = 0.024;
b23 = 0.0014;

b11 = 1 - b12;
b21 = b12*mateq/mueq;
b22 = 1 - b21 - b23;
b32 = b23*mueq/mleq;
b33 = 1 - b32;


positive variables
    bb12          Carbon cycle transition matrix
    bb23
  
    epsMAT1(ttt,n)
    epsMAT2(ttt,n)
    epsMU1(ttt)
    epsMU2(ttt)
    epsML1(ttt)
    epsML2(ttt)
;


VARIABLES
    MAT(ttt,n)          Carbon concentration in atmosphere GtC
    MU(ttt,n)           Carbon concentration in shallow oceans Gtc
    ML(ttt,n)           Carbon concentration in shallow oceans Gtc

    objCO2
;

*****************************************************************

EQUATIONS

    MMAT(ttt,n)          Atmospheric concentration equation
    MMU(ttt,n)           Shallow ocean concentration
    MML(ttt,n)           Shallow ocean concentration

    L1errMAT(ttt,n)
    L1objfunCO2
;

** Equations of the model

MMAT(ttt+1,n)..        MAT(ttt+1,n)    =E= (1-bb12)*MAT(ttt,n) + (MATEQ*bb12/MUEQ)*MU(ttt,n) + Etotal(ttt,n);
MMU(ttt+1,n)..         MU(ttt+1,n)     =E= bb12*MAT(ttt,n) + (1-MATEQ*bb12/MUEQ-bb23)*MU(ttt,n) + bb23*MUEQ/MLEQ * ML(ttt,n);
MML(ttt+1,n)..         ML(ttt+1,n)     =E= bb23*MU(ttt,n) + (1-bb23*MUEQ/MLEQ)*ML(ttt,n);

L1errMAT(ttt,n)..  (MAT(ttt,n)-CA(ttt,n))/CA(ttt,n) =E= epsMAT1(ttt,n)-epsMAT2(ttt,n);
    
L1objfunCO2..     objCO2 =e= SUM((n,ttt), (epsMAT1(ttt,n)+epsMAT2(ttt,n)));
				
**********************************

** Solution options
option iterlim = 99900;
option reslim = 99999;
option solprint = on;
option limrow = 0;
option limcol = 0;
option nlp = conopt;

model L1CalibCO2 / MMAT, MMU, MML, L1errMAT, L1objfunCO2 /;

bb12.l = b12;
bb23.l = b23;
MAT.l(ttt,n)$(ord(ttt)=1) = mat0;
MU.l(ttt,n)$(ord(ttt)=1) = mu0; 
ML.l(ttt,n)$(ord(ttt)=1) = ml0; 
loop(ttt$(ord(ttt)<card(ttt)),
  MAT.l(ttt+1,n) = (1-bb12.l)*MAT.l(ttt,n) + (MATEQ*bb12.l/MUEQ)*MU.l(ttt,n) + Etotal(ttt,n);
  MU.l(ttt+1,n) = bb12.l*MAT.l(ttt,n) + (1-MATEQ*bb12.l/MUEQ-bb23.l)*MU.l(ttt,n) + bb23.l*MUEQ/MLEQ * ML.l(ttt,n);
  ML.l(ttt+1,n) = bb23.l*MU.l(ttt,n) + (1-bb23.l*MUEQ/MLEQ)*ML.l(ttt,n);
);

MAT.fx(ttt,n)$(ord(ttt)=1) = MAT0;
MU.fx(ttt,n)$(ord(ttt)=1) = MU0;
ML.fx(ttt,n)$(ord(ttt)=1) = ML0;

solve L1CalibCO2 minimizing objCO2 using nlp ;

b12 = bb12.l;
b23 = bb23.l;

file CalibPathsCarbon_fitted;
put CalibPathsCarbon_fitted;
CalibPathsCarbon_fitted.nw = 12;
CalibPathsCarbon_fitted.nr = 2;
CalibPathsCarbon_fitted.nz = 1e-12;
loop(n,
loop(ttt,
  put ttt.tl:4:0;
  put CA(ttt,n):14:6;
  put MAT.l(ttt,n):14:6;
  put /;
);
);
